GISデータは 地図で見る統計(統計GIS)のデータダウンロードから入手できる。 ここでは、平成22年度国勢調査(小地域)2010/10/01から男女別人口総数及び世帯総数を選択し、広島県広島市南区のデータを世界測地系緯度経度?Shape形式でダウンロードした。
library(maptools)
pj <- CRS("+proj=longlat +datum=WGS84") # 世界測地系
d <- readShapeSpatial("WGS84/h22ka34103.shp", proj4string = pj) # 町丁?字境界データ
d <- subset(d, d@data$HCODE == 8101) # 水上の境界を除外
plot(d) # 町丁?字境界図
spパッケージのspplotで人口分布の視覚化
library(sp)
(brks <- round(quantile(d@data$JINKO, 0:10/10))) # 分位点で色分け
## 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
## 0 60 287 503 747 908 1075 1270 1592 1912 3526
spplot(d, zcol = "JINKO", col.regions = terrain.colors(length(brks) - 1), at = brks)
GoogleMapに人口で色分けした地図を重ねたhtmファイルを作成しブラウザで表示する。ポリゴン領域をクリックすれば属性情報も表示されるが、日本語は文字化けする。
library(plotGoogleMaps)
(brks <- round(quantile(d@data$JINKO, 0:10/10))) # 分位点で色分け
plotGoogleMaps(d, zcol = "JINKO", layerName = "population", filename = "minamiku.html",
colPalette = terrain.colors(length(brks) - 1), at = brks, strokeColor = 1,
mapTypeId = "ROADMAP")